clear
clc
close all
N=[64 128 256 512 1024];
z1=1;
z2=4;
alp=2;
%h=(z2-z1)/(3*N);
e=3;
ux =@(x) exp(-10.*x.^2);
u1=ux(z1);
syms x;
yu=ux(x);
dyu=diff(yu,x);
gv=yu-alp.*dyu;
fx= str2func(['@(x)',vectorize(string(gv))]);
% [ujn,err,xx,uex]=ode_left_fem(z1,z2,alp,e,N,u1,ux,fx);

for i=1:length(N)
    [ujn,err,xx,uex]=ode_left_fem(z1,z2,alp,e,N(i), u1,ux,fx);
    %error:L2error,L∞error
    h1=(z2-z1)/(N(i));
    
    L2err(i)=norm(err,2)*h1^0.5;
    Lwuerr(i)=norm(err,Inf);
    figure(1); semilogy(xx,err); hold on;
    figure(2);plot(xx,err);
    set(gca,'FontSize',10)
    xlabel('x','FontSize',12)
    ylabel('|u(x_i,1)-u_i^1|','FontSize',12)
    
end
%ORDER
L2order=log(L2err(1:end-1)./L2err(2:end))./log(N(1:end-1)./N(2:end));
Lwuorder=log(Lwuerr(1:end-1)./Lwuerr(2:end))./log(N(1:end-1)./N(2:end));
fprintf('fem,\n')
fprintf('ode\n')
fprintf('u-au_x=f(x) x∈(z_1,z_2)\n')
L_G_output(N,L2err,Lwuerr,L2order,Lwuorder);







function[ujn,err,xx,uex]=ode_left_fem(z1,z2,alp,e,N,u1,ux,fx)
x=linspace(z1,z2,3*N+1);
xx=x';
x1=linspace(z1,z2,N+1);
xx1=x1';
[mm,tt]=basis(e,N,xx1);
m_ha=mm(2:end,:);
t_ha=tt(2:end,:);
min_ha=m_ha(:,2:end);m0_ha=m_ha(:,1);
tin_ha=t_ha(:,2:end);t0_ha=t_ha(:,1);
f_ha=fx(xx);
F=m_ha*f_ha-m0_ha*u1+alp.*(t0_ha*u1);
ujn_in=(min_ha-alp.*tin_ha)\F;
ujn=[u1;ujn_in];
uex=ux(xx);
err=abs(uex-ujn);

end


% m=full(mm);
% t=full(tt);

function [mm,tt] = basis(e,N,xx)
M = [ 8/105,    33/560,  -3/140,   19/1680
      33/560,   27/70,   -27/560,  -3/140
      -3/140,   -27/560, 27/70,    33/560
      19/1680,  -3/140,  33/560,   8/105  ];
  
T = [ -1/2,     57/80,   -3/10,    7/80
      -57/80,   0,       81/80,    -3/10
      3/10,     -81/80,  0,        57/80
      -7/80,    3/10,    -57/80,   1/2    ];
 
e1 = e+1;
% mm = sparse(zeros(e*N+1));  
mm = sparse((e*N+1));  
tt = mm; 

h = xx(2:end) - xx(1:end-1);

for k=0:N-1
    mm(e*k+1:e*k+e1,e*k+1:e*k+e1) = M*h(k+1);
    tt(e*k+1:e*k+e1,e*k+1:e*k+e1) = T;
end
for k=1:N-1
    mm(e*k+1,e*k+1) = M(1,1)*h(k+1) + M(e1,e1)*h(k);
    tt(e*k+1,e*k+1) = T(1,1) + T(e1,e1);
end
mm(1,1)=M(1,1)*h(1);
tt(1,1)=T(1,1);
mm(3*N+1,3*N+1)= M(e1,e1)*h(N);
tt(3*N+1,3*N+1)= T(e1,e1);
end

function  L_G_output(NE,L2err,Lwuerr,L2order,Lwuorder)
fprintf('========================================================================\n')
fprintf(' NE     &L2error      &L∞error     &L2order  &L∞order:\n')
for i=1:length(NE)
    if i==1
       fprintf('%3d     &%5.3e    &%5.3e    &         &\n',NE(i),L2err(i),Lwuerr(i)) 
    else
       formatSpec = '%3d     &%5.3e    &%5.3e    &%5.3f    &%5.3f\n';
fprintf(formatSpec,NE(i),L2err(i),Lwuerr(i),L2order(i-1),Lwuorder(i-1))
%     fprintf('%d     &%5.3f   &%5.3f\n',NE(i),L2order(i),Lwuorder(i))
    end
end
end
